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METHOD OF DETERMINING TEE EQUIVALENT PERMEABILITY OF 
A FRACTURE NETWORK IN A SUB - SURFACE , MULTI - LAYERED 

MEDIUM 

The invention relates to a method of .determining 
the equivalent fracture permeability of a fractured 
network in a sub-surface multi-layered medium, which 
can be used to build a more realistic model of a 
fractured sub- surface geological structure. The method 
can be used by reservoir engineers as a means of 
producing reliable predictions of oil flows, for 
example . 

Fractured reservoirs are an extreme type of 
heterogeneous reservoirs made up of two contrasting 
media, a matrix medium containing the majority of the 
oil and of a low permeability and a fractured medium 
accounting for less than 1% of the oil in place and 
highly conductive. The fractured medium itself may be 
very complex having different fracture sets, each 
characterised by its density, length, orientation, tilt 
and aperture, 3D images of fractured reservoirs can not 
be used directly in the form of input data for 
simulating a reservoir. It has long been considered 
unrealistic to represent fracture networks within 



reservoir flow simulators because the layout of the 
network is partially unknown and because of the 
numerical limitations inherent in juxtaposing numerous 
cells of extremely contrasting sizes and properties. A 
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simplified but realistic way of modelling such media 
would therefore be of great interest to reservoir 
engineers . 

The "dual porosity approach", as taught, for 
example, by Warren J.E. et al in "The Behaviour of 
Naturally Fractured Reservoirs M , SPE Journal (September 
1963),- 245-255, is known in the art as a means of 
interpreting single -phase flow behaviour observed when 
testing a fractured reservoir. In accordance with this 
basic model, any elementary volume of the fractured 
reservoir is modelled in the form of a set of identical 
parallelepipedic blocks limited by an orthogonal system 
of continuous uniform fractures oriented in the 
direction of one of the three main flow directions. On 
the reservoir scale, fluid flows through the fractured 
medium only and fluid exchanges occur locally between 
the fissures and the matrix blocks. 

Numerous simulators have been developed for 
fractured reservoirs using a model of this type with 
specific improvements in respect of the way flows are 
exchanged between matrix and fracture, governed by 
capillary, gravitational and viscous forces and 
compositional mechanisms as well as matrix- to-matrix 
flow exchanges (dual -permeability dual-porosity 
simulators) . Numerous examples of prior art techniques 
are cited in the reference works listed below. 
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Thomas, Ij.K. et al : "Fractured Reservoir 
Simulation", SPE Journal (February 1983) 52-54; 

Cuandalle, P. et al : "Typical Features of a New 
Multipurpose Reservoir Simulator", SPE 16007, presented 
at the 9th Symposium on reservoir simulation held in 
San Antonio, Texas,, 1-4 February 1987; 

Coats, K. H. : "Implicit Compositional Simulation 
of Single -Porosity and Dual-Porosity Reservoirs", SPE 
18427 presented to the SPE Symposium on reservoir 
simulation held in Houston, Texas, 6-8 February 1989. 

One of the problems which reservoir engineers 
encounter is that of applying parameters to this basic 
model in order to produce reliable flow predictions. 
The basic petro-physical properties of fractures and 
matrices and the size of the matrix blocks in 
particular must be ascertained for each cell of the 
flow simulator- Whilst the permeability of the matrix 
can be estimated from cores, there is no simple method 
of estimating the permeability of the fracture network 
contained within the cell, i.e. the equivalent fracture 
permeability and it is necessary to take account of the 
geometry and properties of the fracture network. 

One direct method is known whereby steady- state 
flow is determined in a fracture network. It consists 
in using conventional fine and regular grids to divide 
up both the fractures and the matrix blocks of the 
parallelepipedic rock volume considered. For various 
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reasons, this known method does not produce reliable 
results unless the volume of fractured rock is divided 
using a grid with an extremely high number of cells, 
which requires enormous computing resources . 
5 Other specific models which compute equivalent 

permeabilities of 2D or 3D fracture networks, for 
example, are also known from: 

Odling, N.B. : "Permeability of Natural and 
Simulated Fracture Patterns" , Structural and Tectonic 
10 Modelling and its Application to Petroleum Geology, NPF 
Special Publication 1, 365-360, Elsevier, Norwegian 
Petroleum Society (NPF) 1992; 

- Long, J.C.S. : et al : "A Model for Steady Fluid 
Flow in Random Three -Dimensional Networks of Disc- 

15 Shaped Fractures", Water Resources Research {August 
1985), vol. 21, 8, 1105-1115; 

- Cacas, M.C. et al : "Modelling Fracture Flow with . *.-.r- 
a Stochastic Discrete Fracture Network r Calibration . 
and Validation. 1. The Flow Model 1 * , Water Resources 

20 Research CMarch 1990) vol. 26, 3 ; 

Billaux, D. : "Hydrogeologie des milieux 
fractures. Geometrie, connectivity et comportement 
hydraulique" , Doctoral thesis presented to the Ecole 
Nationale Superieure des Mines de Paris, BRGM document 
25 186, Editions du BRGM, 1990; 

Robinson, P.C.: "Connectivity, Flow and 
Transport in networks, Models of Fractured Media", 
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Doctoral thesis, St Catherine's College, Oxford 
University, Ref . : TP1072, May 1984. 

The present invention relates to a method of 
determining the equivalent fracture permeability of a 
5 fracture network in a sub- surf ace, multi- layered 
medium . 

The method is characterised in that it consists of 
the following steps: 

- the fracture network is split up into fracture 
10 elements (such as rectangles, for example) and nodes 

are defined to represent inter-connected fracture 
elements in each layer of the medium and 

- fluid flows through the divided network are 
determined by imposing boundary pressure conditions and 

15 fluid transmissivities to each pair of adjacent nodes. 

More specifically, the method is characterised in 

that: 

- the medium is divided into a series of parallel 
layers, each extending in a reference plane 

20 perpendicular to the reference axis and defined by a 
coordinate on the said axis, 

each fracture is divided into a series of 
rectangles limited by two adjacent layers along the 
said reference axis and the rectangles are itemised by 

25 attributing geometric and physical characteristics to 
them such as coordinates and sizes of the rectangles 
and hydraulic conductivities of the fractures. 



snt By,: 




18006661233; 



Mar 



00 9:41 ; . 



Page 14/22 



6 



10 



15 



20 



- nodes are placed in each layer for all the 
inter- connected fractures, and 

transmissivity factors are computed and flow 
equations solved for all the pairs of adjacent nodes in 
order determine the equivalent permeabilities of the 
medium in three orthogonal directions. 

In a preferred embodiment, the equivalent 
permeability of the medium includes directly 
determining the equivalent permeability anisotropy 
tensor and calibrating absolute permeability values 
using well test results. 

The method as summarised allows characteristic 
models of fractured reservoirs to be systematically 
linked with dual porosity simulators in order to 
produce a more realistic model of a sub-surface, 
fractured geological structure. The method can be ' 
implemented by reservoir engineers in the field of ~- 
petroleum production as a means of providing reliable • ■ ■■ 
flow predictions . 

Other features and advantages of the method of the 
invention will become clearer from the following 
description of embodiments, given by way of example and 
not restrictive in any respect, and with reference to 
the appended drawings, in which: 

figure l shows a 3D image of a fracture network, 
for example, stochastically generated from observations 
of and measurements taken on a sandstone outcrop, 
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figures 2, 3 illustrate a fracture divided into a 
series of rectangles R, 

figure 4 shows the structure of input data 
specifying the fracture attributes , 

figure 5 shows the preferred method of dividing up 
a fracture plane, 

figures 6, 7, 8 and 9 provide a schematic 
illustration of how transmissivity factors are computed 
for different positions of nodes relative to one 
another or with a boundary. 

The equivalent permeabilities of a 3D fracture 
network are determined below using a numerical 
technique based on the known "resistor network" method, 
described for example in the publication by Odling, 
N.E., mentioned above. With this method, the matrix is 
assumed to be impermeable so as to be consistent with 
the dual porosity approach. In reservoir simulators, 
matrix-to-fracture and matrix- to -matrix flows are 
actually computed separately from flows within the 
fractures. 

The 3D fracture network considered is assumed to 
represent, in a volume equal to a reservoir cell, the 
real distribution of fractures produced by integrating 
fracture attributes of the field in a characterisation 
model. The main purpose of single-phase flow 
computations in the 3D fracture network is to estimate 
the equivalent permeability anisotropy (Kv/Kh and 
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Ky/Kx) of the fracture cell considered, which is an 
important parameter governing the behaviour of multi- 
phase flows in reservoirs. The equivalent permeability 
values obtained by these computations could in practice 
be compared with the results of well test as a means of 
calibrating fracture attributes such as fracture 
hydraulic conductivities {or equivalent hydraulic 
apertures) , which a priori can be defined to a mediocre 
standard only. 

The equivalent permeability results can also be 
used to determine a permeability tensor, the main 
directions of which will enable optimum orientation of 
the reservoir model grid. However, specific boundary 
conditions are needed in order to obtain this 
information. Lateral no- flow boundary conditions 
imposed on the four lateral faces of the 
parallelepipedic volume under study do not give access 
to non- diagonal terms of the equivalent permeability 
tensor whereas linearly varying potentials (or. 
pressures) across the lateral faces will allow the 
direction of potential gradient within the anisotropic 
medium to be imposed and the non-diagonal permeability 
terms of lateral flows to be derived directly. 

The techniques used to integrate natural 
fracturing data into fractured reservoir models are 
well known. Fracturing data are essentially geometric 
in nature and incorporate measurements of the density, 




18006661233; 



Mar 




00 9:42; 



. 9 



length, azimuth and tilt of fracture planes either as 
observed from outcrops, mine drifts or cores or 
inferred from well logging. Different fracture sets can 
be distinguished and characterised by different 
statistical distributions of the fracture attributes. 
Once the fracture patterns have been characterised, 
numerical networks of these fracture sets can be 
created using a stochastic method consistent with the 
statistical distributions of the fracture parameters . 
Such methods are described in patents FR-A-2 . 725 . 814, 
2.725.794 or 2.733.073 filed by the applicant, for 
example. 

The method of the invention is applied to images 
of fractured geological structures of various sizes or 
volumes and/or at various locations which are generated 
by a fracture model. Figure 1 illustrates such an 
image. . 

INPUT DATA 

Before developing the procedures recommended as a 
means of determining the equivalent hydraulic 
parameters of 3D fracture images, a major step consists 
in defining a common structure for the input data to be 
applied to these images so that they can be processed 
independently of the processing tool used to generate 
them. 

As shown in figures 2, 3, it is assumed that the 
fractures F are essentially vertical (i.e. 
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perpendicular to the layer limits) . However, a same 
data structure can be applied to fractures which 
deviate slightly from the vertical. The 3D image is 
divided vertically, complying with the , actual 
5 geological layering if this information is available. 
Otherwise, an arbitrary division is applied to the 
image.- Each horizontal layer L is characterised by its 
vertical coordinate zl» in the reference system of 
coordinates (OX, OY, 02) . 

10 A series of rectangles "R must be defined for each 

layer L. Each rectangle consists of a fracture plane 
element between the limits of a given layer. 
Consequently, each natural fracture consists of a 
series of superposed rectangles R and is assigned an 

15 origin (origin of the fracture) . Each rectangle is 
defined by: 

- the three coordinates {xO, yO, zO) of the origin 
of the rectangle O. All the origin points of the 
rectangles making up any given natural fracture are 

20 located on the same vertical line (or highest dip)* 
plotted from the origin of the fracture; 

- the coordinates of the horizontal unit, vector 
. (xH, yH) and the vertical unit vector (xV, yV) 

defining the orientation of the rectangle within the 
25 reference system of coordinates, where xVertical and 
yVertical are equal to zero in the case of vertical 
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fractures but are considered as input data in order to 
be able to deal with non- vertical fractures ,- 

- the two algebraic horizontal lengths 1- and 1+ 
separating the origin of the rectangle and the two 
lateral (vertical} limits of this rectangle ,- 

- the height h of the rectangle, i.e. the length 

— ► 

of the rectangle in the direction . , which is the 
thickness of the layer if the division in direction 



Darcy's law on fracture flow (for a pressure gradient 



— , the flow rate in the fracture of a height h is 



where k » a 2 /12 (using Poiseuille's idealised 
representation of fractures) is the intrinsic 
permeability of the fracture and a its equivalent 
hydraulic aperture. The hydraulic conductivity c is a 
reference value given for a direction of the maximum 
horizontal stress parallel with the direction of the 
fracture ; 

- the two upper and lower adjacent rectangles UR, 



corresponds to the geology; 



the hydraulic conductivity c derived by applying 




LR; 



the fracture series FS to which the considered 



rectangle belongs; 
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- the angle of orientation ct e of the direction of 
maximum horizontal stress taken from the axis (OX) 
within the reference system of coordinates ,- 

- for each series of fractures, a correlation 
5 table correlating 1) the angle between the direction of 

maximum horizontal stress and the direction of the 
fracture (azimuth) with 2) the hydraulic conductivity c 
of the equivalent hydraulic aperture a previously 
defined. The terms "horizontal M and "vertical" used in 

10 this context relate to the respective directions 
parallel with and perpendicular to the layer limits, 
which are assumed to be horizontal in this case . The 
layer limits divide the fracture plane in the 
"vertical" direction. It should be pointed out that the 

15 above-mentioned input data 1) are suitable for all the 
existing software tools used to characterise and 
generate fractures and 2) could be used to divide up a 
network of fractures slightly off from the vertical,, 
i.e. which are not perpendicular to the layer limits. 

20 OPERATING PROCEDURES 

The operating procedures and validation tests for 
the method used to compute the permeability anisotropy 
of the fracture network taken as a whole are described 
below using the 3D image coded in this way. The 

25 numerical procedure used to compute the equivalent 
permeabilities of a 3D fracture network will also be 
described. 
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The problem consists in finding the distribution 
of flow rates in the network for the following boundary 
conditions imposed on the limits of the 
parallelepipedic volume studied, i.e. fixed pressures 
5 imposed on two opposing faces and pressures varying 
linearly across the four lateral faces (between the 
values imposed on the other two faces) . 
The main steps are set out below: 
1) Dividing up the network 

10 Using the definitions given for the input data 

structure, the fracture network is divided up into a 
series of "nodes" N, each node being positioned at the 
centre of the intersecting segments IS of two 
rectangles R <i,e. of two fracture planes within a 

15 given layer) . As illustrated in figure 5, additional 
nodes AN are positioned above and below the preceding 
nodes N in order to represent other rectangles dividing 
up the fractures and minimising the flow lengths within 
a given fracture. In figure 5, BL is a lateral limit of 

20 two adjacent fracture cells. 

Once the network has been divided up, a screening 
procedure is applied to this fracture network in order 
to eliminate isolated nodes or groups of nodes with no 
link to one of the lateral limits FL of the 3D volume 

25 studied since the fissures "screened" in this way do 
not contribute to fluid transport and may impede the 
computational procedures used to find the pressures at 
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fracture nodes during a steady- state flow through the 
network . 

2) Calculating transmissivities 

A transraissivity factor T is calculated for each 
pair of connected nodes using the relationship: 



where c is the hydraulic conductivity of the 
fracture, k the intrinsic permeability of the fracture, 
a the aperture of the fracture, h the height of the 
fracture and 1 the distance between two nodes of a 
fracture . 

Different situations have to be considered 
depending on the respective position of the two nodes. 
For nodes contained within the same layer (fig. 6) the 
horizontal transmissivity factor T is obtained directly 
in the form of the distance (11+12) separating the two 
nodes in the flow direction (fig. 7) . in the case of 
nodes located in two different layers (fig. 9) , the 
horizontal transmissivity factor is the arithmetic sum 
of the transmissivity factors (T 1 +T") relating to the 
two fracture plane elements of the superposed fracture 
cells. It comprises a flow length equal to half the 
sum of the two layer thicknesses hi and h2. For the 
additional nodes as previously defined, a single 
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transmissivity factor is calculated for this fracture 
p 1 ane el emen t . 

The transmissivity factor T between a node and a 
limit of the 3D volume studied is expressed in the same 
way as between two nodes, depending on one of the 
following two situations. 

In the case of a lateral vertical limit, the 
transmissivity factor T can be expressed directly for a 
single fracture plane element {fig. 8) and takes the 
form of the sum of two transmissivities if the node and 
the limit are linked by two fracture planes. 

In the case of a horizontal upper or lower limit, 
the vertical transmissivity factor can be expressed by 
taking a flow length equal to half the thickness of the 
layer (fig. 9) . 

3) Flow equations 

In steady state, an incompressible single -phase 
flow through the fracture network is determined by 
solving a series of n equations, one for each node, 
well known in the art. Each equation expresses the fact 
that the total flow rate is zero at each fracture node. 
In order to calculate a permeability tensor, a constant 
pressure is imposed on each of the limits upstream and 
downstream. A pressure is imposed which varies linearly 
as a function of the position between the upstream and 
downstream limits. 



ent By: - ; 18006661233; Mar^6-00 9:44; Page 2 

16 



The equivalent permeability matrix (Ki j ) 
determined previously is diagonalised in order to 
calculate the main flow directions with the respective 
equivalent permeabilities in these directions. 
5 In practice, the problem is often limited to 

finding the main horizontal flow directions U and V 
given that the direction perpendicular to the layer 
limits (generally vertical) is always taken as the z 
axis. In such a case, only the extra -diagonal terms 
10 and need to be calculated: they can be obtained 

with the following mixed boundary conditions : 

- horizontal flows are computed with impermeable 
top and bottom faces and pressures varying linearly on 
the vertical faces parallel with the flow direction; 
15 - vertical flow is computed with all lateral faces 

being impermeable. 

A simplified permeability tensor is thus obtained, 
from which the main horizontal flow directions CJ and V 
can easily be derived: 



20 



^00 K^J 



Validation 

The method has been successfully validated against 
25 the above-mentioned reference single -phase flow 
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computations performed with a conventional reservoir 
simulator. The reference computations were produced on 
fine regular grids dividing up the fractures as well as 
the matrix blocks of the parallelepipedic fractured 
5 rock volume considered. For a given flow direction, a 
fixed injection pressure and production pressure were 
imposed on the inlet and outlet faces and the resulting 
flow rate was calculated for conditions under which 
there is no lateral flow. 
10 Three steps were performed to validate the 

computation of: 

- the equivalent vertical permeability of a volume 
of rock crossed by a single fracture, this latter being 
represented by several nodes corresponding to 

15 intersections with small disconnected fractures; 

- the equivalent horizontal permeabilities (in a 
2D flow geometry) and the main flow directions; 

- the equivalent permeabilities and permeability 
anisotropy in a single network involving 3D flow 

20 geometry. 

The results obtained for the third step {for a 3D 
flow geometry) are set out in the table below. A 
reference analytical solution can also be computed for 
the horizontal flow directions since the flow geometry 

25 is a 2D flow geometry in these directions (the 3D flow 
geometry relates to the z direction) , 
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Equivalent 
permeabilities 
(md) 


FINE GRID 
simulation 


PRESENT 
METHOD 


ANALYTICAL 
solution 


Kx 


0.119 


0.120 


0.120 


Ky 


0.224 


0.227 


0.226 


K2 


0.255 


0.267 




An i sot ropy 
Kz/(KxKy) DS 


1.56 


1.62 





The results produced by the present method are 
clearly very close to the corresponding values obtained 
using the analytical solution and the fine grid 
5 solution for the directions X and Y. 

Furthermore, the difference between the vertical 
equivalent permeability values relating to 3D flow is 
still acceptable. Hence, the anisotropy ratio, which is 
1.6, is satisfactorily predicted by the method with a 
10 very limited number of cells. 

The method of the invention, which provides a 
readily transposed representation of a natural fracture 
network is well suited to fracture flow computations. 
It may also prove useful in improving the original 
15 image of the fracture network. Such an image is in fact 
obtained from a stochastic fracture generator using as 
input data the results produced by integrating filed 
fracture data in a fracture characterisation model as 
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described in patents FR-A-2.72S.S14, 2.725.794 or 
2.733.073 filed by the applicant, mentioned above. Once 
divided up by the method of the invention, such image* 
can easily be modified to fit the geological rules . For 
example, the systematic interruption of a given 
fracture against any other series of fractures can be 
taken into account in the original image by eliminating 
fracture plane elements from a given set extending 
beyond the intersected fractures of the other set. 
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CULIMS 

1. A method for determining the equivalent 
fracture permeability of a fracture network in a sub- 
surface, fractured, multi-layered medium from a 
predetermined representation of said network, 
comprising the steps of: 

dividing each fracture <F) in the fracture 
network into fracture elements and defining ■ nodes (N> 
representing inter -connected fracture elements in each 
layer of the medium and 

determining fluid flows through the fracture 
network, imposing boundary pressure conditions, and 
fluid transmissivities for each pair of adjacent nodes. 

2. A method as claimed in claim 1, characterised 
in that: 

the medium is divided into a series of 
parallel layers each extending ^in a reference plane 
perpendicular to the reference axis and each defined by 
a coordinate on said axis, 

each fracture is divided into a series of 
rectangles limited on said reference axis by two 
adjacent layers and the rectangles are itemised by 
attributing geometric and physical characteristics to 
them such as rectangle coordinates and sizes and 
hydraulic conductivities of the fractures, 



ent By: - ; 18006661233; Mar^6-00 9:45; Page 3 



21 



nodes are placed in each layer for all che 
inter- connected fractures and 

transmissivity factors are computed and flow 
equations solved for all the pairs of adjacent nodes in 
5 order to determine the equivalent permeabilities of the 
medium in three orthogonal directions. 

3. A method as claimed in one of the preceding 
claims, characterised in that the equivalent 
permeability of the medium includes directly 

10 determining an equivalent permeability tensor and 
calibrating absolute permeability values from well test 
results. 

4 . A method as claimed in one of the preceding 
claims, characterised in that said geometric and 

15 physical attributes are the coordinates and dimensions 
of the rectangles and the hydraulic conductivities of 
the fractures. 

5. A method as claimed in any preceding claim 
and as hereinbefore described with reference to the 

2 0 accompanying drawings. 
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